from copy import copy
import re
import math
from number_theory import multinomial_coefficients, bitcount, newton_method_sizes
rpm = re.compile('[+-]')


class RPolyOverflowError(OverflowError):
  pass

class BaseRPoly:
  """
  abstract base class for polynomial rings of a ring K
  subclasses:
  RPoly polynomial ring on a commutative ring
        with no elements with more than one term
  MRPoly polynomial ring on a commutative ring
        with elements with more than one term
  NCPoly polynomial ring with noncommutative base ring

  The objects rp in one of these polynomial rings
  construct elements of Poly, a multivariate polynomial ring with
  monomial ordering.

  rp.field generates the elements in K
  K(0)  zero element in K
  K(1)  unit element in K
  
  Order can be lex, grlex, or a tuple indicating the weight of the
  degrees (weight order, see [CLO] D. Cox, J. Little, D.O'Shea 
  'Ideals, varieties and algorithms' p74), 
  with lex order breaking ties.
  The exponents of a monomial are packed in a python integer.
  There is an extra bit for each variable, to allow checking 
  overflow in multiplication and to check divisibility.
  expv1, expv2 encoded exponents of monomials m1, m2
  m3 = m1*m2 ; if expv3 = expv1 + expv2 satisfies expv3&maskof == 0,
  then expv3 encodes the exponents of m3.
  m1 divides m2 iff (expv2 - expv1)&maskof == 0; expv2 - expv1 encodes
  the exponents of m2/m1
  For the grlex and weight orders there is an extra hidden variable, 
  to keep track of the weight degree of the monomial, which is the
  total degree in the case of grlex.
  In pol_gens the variables are ordered in increasing order.
  
  The string representation of polynomials is in decreasing order;
  in it ^ instead of ** is used as exponent symbol.

  Note: in the doctests Fraction is used; using gmpy.mpq is faster.
  Replacing
  from fractions import Fraction as mpq
  with
  from gmpy import mpq
  the output is the same
  >>> from rmpoly import *
  >>> from fractions import Fraction as mpq
  >>> rp = RPoly(['y','x'],6,mpq)
  >>> y,x = rp.gens()
  >>> x + y**2
   +x +y^2
  >>> rp = RPoly(['y','x'],6,mpq,order='grlex')
  >>> y,x = rp.gens()
  >>> x + y**2
   +y^2 +x
  >>> rp = RPoly(['y','x'],6,mpq,order=(3,1))
  >>> y,x = rp.gens()
  >>> x + y + x**2 + x**3
   +x^3 +y +x^2 +x
  """
  def __init__(self, pol_gens, bits_exp, field, **kwds):
    if isinstance(pol_gens,basestring):
      pol_gens = pol_gens.split(',')
    self.pol_gens = tuple(pol_gens)
    self.ngens = len(pol_gens)
    self.bits_exp = bits_exp
    self._bits_exp = int(bits_exp+1)
    self.field = field
    if 'base_field' in kwds:
      self.base_field = kwds['base_field']
    else:
      self.base_field = None
    if 'order' in kwds:
      self.order = kwds['order']
    else:
      self.order = 'lex'
    order = self.order 
    self.mask_exp = (1<<self._bits_exp)-1
    if order == 'grlex' or isinstance(order, tuple):
      self.ordertype = 1
      self.mask_vars = (1<<((self.ngens)*self._bits_exp)) - 1
      if order == 'grlex':
        self.ordertuple = tuple([1]*self.ngens)
      else:
        self.ordertuple = map(int,order)
    else:
      # order='lex'
      if order:
        assert order == 'lex'
      self.ordertype = 0
    maskof = 0
    mask1 = (1<<(self._bits_exp-1))
    for i in range(self.ngens):
      maskof += mask1<<(i*self._bits_exp)
    self.maskof = maskof
    self.gens_dict = {}
    for i in range(self.ngens):
      self.gens_dict[self.pol_gens[i]] = i
    if self.ordertype == 1:
      self.gens_dict['_t'] = self.ngens
    if 'parens' in kwds:
      self.parens = kwds['parens']
    else:
      self.parens = False
    if 'commuting' in kwds:
      self.commuting = kwds['commuting']
    else:
      self.commuting = True
    if isinstance(field,BaseRPoly):
      self.parens = True
      if not field.commuting:
        self.commuting = False
    str_field = str(self.field)
    if 'mpc' in str_field or 'Complex' in str_field:
      self.parens = True

  def __str__(self):
    try:
      field_name = self.field.__name__
    except:
      field_name = str(self.field)
    s = 'RPoly with ngens=%d bits_exp=%d field=%s' % \
        (self.ngens, self.bits_exp, field_name)
    return s

  def gens(self):
    """return the list of the variables
    >>> from rmpoly import *
    >>> from fractions import Fraction
    >>> rp,y,x = rgens('y,x',6,Fraction)
    >>> (x+y)**2
     +x^2 +2*x*y +y^2
    """
    bits_exp = self._bits_exp
    ordertype = self.ordertype
    if ordertype:
      ordertuple = self.ordertuple
      shift = self.ngens*bits_exp
    a = []
    one = self.field(1)
    for i in range(self.ngens):
      p = Poly(self)
      expv = 1<<(i*bits_exp)
      if ordertype:
        expv += ordertuple[i]<<shift
      p[expv] = one
      a.append(p)
    return a

  def code_expv(self, a):
    """convert the tuple a of exponents for monomial x0^a[0]*x1^a[1]*...
       in encoded expv; bits are added to expv in grlex and weight
       orders to keep track of the
       the weight degree of the monomial sum(a[i]*ordertuple[i])
    """
    bits_exp = self._bits_exp
    expv = 0
    if self.ordertype:
      deg = 0
      ordertuple = self.ordertuple
      for i, n in enumerate(a):
        deg += n*ordertuple[i]
        expv += n<<(i*bits_exp)
      expv += (deg<<(self.ngens*bits_exp))
    else:
      for i, n in enumerate(a):
        expv += n<<(i*bits_exp)
    return expv

  def decode_expv(self, expv):
    """decode expv to [e0,e1,..] exponents of a monomial"""
    bits_exp = self._bits_exp
    mask_exp = self.mask_exp
    a = []
    i = 0
    # throw away the bits for the hidden variable
    if self.ordertype == 1:
      expv &= self.mask_vars
    while(expv):
      exp = expv&mask_exp
      expv >>= bits_exp
      a.append(exp)
    return a

  def new(self,p):
    """
    return a polynomial with RPoly self, with new order or new bits_exp
    >>> from rmpoly import *
    >>> from fractions import Fraction as mpq
    >>> rp,y,x = rgens('y,x',6,mpq)
    >>> p = x + y**2; p
     +x +y^2
    >>> rp1 = RPoly(['y','x'],10,mpq,order='grlex')
    >>> rp1.new(p)
     +y^2 +x
    """
    p1 = self('0')
    decode_expv = p.rp.decode_expv
    code_expv = self.code_expv
    for expv in p:
      degrees = decode_expv(expv)
      expv1 = code_expv(degrees)
      p1[expv1] = p[expv]
    return p1

  def read_monom(self,s):
    """compute expv for the string 'x_0^e_0*...'
    where expv = rp.code_expv((e0,...))
    """
    mask_exp = self.mask_exp
    bits_exp = self._bits_exp
    gens_dict = self.gens_dict
    if self.ordertype:
      ordertuple = self.ordertuple
    field = self.field
    s = s.strip()
    a = s.split('*')
    expv = 0
    # e.g. a = ['x0^6', 'x1^6', 'x10^2']
    # split each element in ind,pw
    if self.ordertype == 1:
      deg = 0
      for x in a:
        t = x.split('^')
        ind = gens_dict[t[0]]
        if len(t) == 2:
          pw = int(t[1])
        else:
          pw = 1
        expv += pw<<(ind*bits_exp)
        deg += pw*ordertuple[ind]
      expv += (deg<<(self.ngens*bits_exp))
      return (expv, coeff)
    else:
      for x in a:
        t = x.split('^')
        ind = gens_dict[t[0]]
        if len(t) == 2:
          pw = int(t[1])
        else:
          pw = 1
        expv += pw<<(ind*bits_exp)
    return expv

  def mon_eval(self,s):
    """compute the tuple (expv,c) for a string 'c*x_0^e_0*...'
    where expv = rp.code_expv((e0,...))
    """
    mask_exp = self.mask_exp
    bits_exp = self._bits_exp
    gens_dict = self.gens_dict
    if self.ordertype:
      ordertuple = self.ordertuple
    field = self.field
    s = s.strip()
    if s[0].isdigit() and '*' not in s:
      return (0,field(s))
    if s[0].isdigit():
      a = s.split('*',1)
      coeff = field(a[0])
      a = a[1].split('*')
    else:
      a = s.split('*')
      coeff = field(1)
    expv = 0
    # e.g. a = ['x0^6', 'x1^6', 'x10^2']
    # split each element in ind,pw
    if self.ordertype == 1:
      deg = 0
      for x in a:
        t = x.split('^')
        ind = gens_dict[t[0]]
        if len(t) == 2:
          pw = int(t[1])
        else:
          pw = 1
        expv += pw<<(ind*bits_exp)
        deg += pw*ordertuple[ind]
      expv += (deg<<(self.ngens*bits_exp))
      return (expv, coeff)
    else:
      for x in a:
        t = x.split('^')
        ind = gens_dict[t[0]]
        if len(t) == 2:
          pw = int(t[1])
        else:
          pw = 1
        expv += pw<<(ind*bits_exp)
      return (expv, coeff)

  def read(self, s):
    """
    return the polynomial of s, in case of polynomial over polynomial
    s string of the form +(coeff1)*m1 + ... +(coeff0)
    where coeff0,coeff1,.. are polynomials in the field self.field
    m1,..  monomial formed by the variables rp.gens()
    The form of the string must be exactly as above; each coeff is in brackets,
    there is always a + in front of the left bracket, if there is a monomial
    multiplying the monomial the right bracket is followed by * (no space allowed);
    the only exception is when s is a number.
    """
    s = s.strip()
    p = Poly(self)
    rp1 = self.field
    if '(' not in s:
      p = Poly(p.rp)
      cp = rp1(s)
      if cp:
        p[0] = cp
      return p
    assert s[0] == '+' and s[1] == '('
    # nparens = number of '(' parenthesis - number of ')' parenthesis
    nparens = 1
    # prev_pos starting position of the subexpression; the subexpression ends
    # when nparens becomes zero.
    prev_pos=2
    ns = len(s)
    pos = 2
    while 1:
      if s[pos] == '(':
        nparens += 1
      elif s[pos] == ')':
        nparens -= 1
      if nparens == 0:
        # found subexpression
        s1 = s[prev_pos:pos]
        # evaluate the subexpression
        cp = rp1(s1)
        if cp:
          pos += 1
          # if s1 is followed by *, then the monomial multiplying it follows
          if pos < ns and s[pos] == '*':
            pos += 1
            # pos1= position of next coefficient subexpression
            pos1 = s.find('+(',pos)
            # +(
            prev_pos = pos1+2
            if pos1 == -1:
              expv = self.read_monom(s[pos:])
            else:
              expv = self.read_monom(s[pos:pos1])
          else:
            expv = 0
            pos1 = -1
          p[expv] = cp
          if pos1 == -1 or pos >= ns:
            break
          else:
            pos = pos1 + 1
          if pos >= ns:
            break
        else:
          # the coefficient is zero, go to the next term
          if pos < ns:
            pos1 = s.find('+(',pos)
            if pos1 == -1:
              break
            prev_pos = pos1+2
            pos = prev_pos
      else:
        pos += 1
        if pos >= ns:
          break

    if nparens != 0:
      raise ValueError, 'nparens=%d' % nparens
    return p

  def __call__(self, s):
    """
    generate a polynomial from a string; the string must be in
    the form c_0*m_0 + c_1*m_1 + ...
    where c_i are coefficients and m_i are monomials
    >>> from rmpoly import *
    >>> from fractions import Fraction as mpq
    >>> rp,y,x = rgens('y,x',6,mpq)
    >>> rp('1/2 + 2/3*x^4 + 1/5*x*y^2')
     +2/3*x^4 +1/5*x*y^2 +1/2
    >>> rp,x = mrgens('x',6,complex)
    >>> p = 1 + 1j + x*(1-1j)
    >>> p1 = rp(str(p))
    >>> p1 == p
    True
    """
    if isinstance(s, Poly):
      if s.rp == self:
        return s
      else:
        raise ValueError('cannot convert polynomial')
    if not isinstance(s, str):
      try:
        c = self.field(s)
        p = Poly(self)
        if c:
          p[0] = c
        return p
      except:
        raise ValueError('s cannot be converted to self.field')
    if self.parens:
      return self.read(s)
    s = s.strip()
    a = []
    if (s[0] == '+' or s[0] == '-'):
      sgn = s[0]
    else:
      sgn = '+'
      s = '+' + s
    it = rpm.finditer(s)
    it.next()
    s1 = s[1:]
    if rpm.search(s1):
      start = 1
      for match in it:
        t=match.span()
        a.append((sgn,s[start:t[0]]))
        start=t[0]+1
        sgn = s[t[0]]
      a.append((sgn, s[t[1]:]))
    else:
      a.append((sgn,s1))
    p = Poly(self)
    for sgn, s1 in a:
      m = self.mon_eval(s1)
      if sgn == '-':
        m = (m[0],-m[1])
      if m[0] in p:
        p[m[0]] += m[1]
      else:
        p[m[0]] = m[1]
    p.strip()
    return p

  def from_mon(self, a):
    """polynomial from the monomial coeff*x_i^j
       a = (i,j,coeff)
    """
    bits_exp = self._bits_exp
    p = Poly(self)
    i = a[0]
    if i.__class__ == str:
      i = self.gens_dict[i]
    j = a[1]
    c = a[2]
    nb = i*bits_exp
    expv = j<<nb
    if self.ordertype:
      expv += (j<<(self.ngens*bits_exp))*self.ordertuple[i]
    p[expv] = c
    return p

  def from_expv(self,expv):
    "monomial from encoded exponent expv"
    p = Poly(self)
    p[expv] = 1
    return p

  def lcm_expv(self,expv1,expv2):
    """ expv1, expv2 encoded exponents for monomials m1,m2
    output expv for lcm(m1,m2)
    """
    mask_exp = self.mask_exp
    bits_exp = self._bits_exp
    expv = 0
    if self.ordertype:
      ordertuple = self.ordertuple
      deg = 0
      for i in range(self.ngens):
        mask = mask_exp<<(i*bits_exp)
        m = max(expv1&mask,expv2&mask)
        deg += (m>>(i*bits_exp))*ordertuple[i]
        expv += m
      expv += deg<<(self.ngens*bits_exp)
    else:
      for i in range(self.ngens):
        mask = mask_exp<<(i*bits_exp)
        m = max(expv1&mask,expv2&mask)
        expv += m
    return expv

class RPoly(BaseRPoly):
  def __init__(self, pol_gens, bits_exp, field, **kwds):
    BaseRPoly.__init__(self, pol_gens, bits_exp, field, **kwds)
    self.__name__ = 'RPoly'


class MRPoly(BaseRPoly):
  """class of polynomials on a ring with elements with more than one term

  >>> from rmpoly import *
  >>> from fractions import Fraction as mpq
  >>> rp1,x = rgens('x', 8, mpq)
  >>> rp2,y = rgens('y', 8, rp1)
  >>> p = x + y + y*x; p
   +( +x +1)*y +( +x)
   >>> isinstance(rp2,MRPoly)
   True
  """

  def __init__(self, pol_gens, bits_exp, field, **kwds):
    BaseRPoly.__init__(self, pol_gens, bits_exp, field, **kwds)
    self.parens = True
    self.commuting = True
    self.__name__ = 'MRPoly'  
  

class NCRPoly(BaseRPoly):
  """class of polynomials on noncommuting ring
  """
  def __init__(self, pol_gens, bits_exp, field, **kwds):
    BaseRPoly.__init__(self, pol_gens, bits_exp, field, **kwds)
    self.parens = True
    self.commuting = False
    self.__name__ = 'NCRPoly'  
  

def rgens(pol_gens, bits_exp, field, **kwds):
  """
  factory function to generate RPoly object and its generators

  >>> from rmpoly import *
  >>> from fractions import Fraction as mpq
  >>> rp,y,x = rgens('y,x',6,mpq)
  >>> x + y**2
   +x +y^2
  >>> rp,y,x = rgens('y,x',6,mpq,order='grlex')
  >>> x + y**2
   +y^2 +x
  >>> rp,y,x = rgens('y,x',6,mpq,order=(3,1))
  >>> x + y + x**2 + x**3
   +x^3 +y +x^2 +x
  """
  rp = BaseRPoly(pol_gens, bits_exp, field, **kwds)
  if rp.parens:
    rp.__class__ = MRPoly
    rp.__name__ = 'MRPoly'
  else:
    rp.__class__ = RPoly
    rp.__name__ = 'RPoly'
  return [rp] + rp.gens()

def mrgens(pol_gens, bits_exp, field, **kwds):
  """
  factory function to generate MRPoly object and its generators
  """
  rp = MRPoly(pol_gens, bits_exp, field, **kwds)
  return [rp] + rp.gens()

def ncrgens(pol_gens, bits_exp, field, **kwds):
  """
  factory function to generate NCRPoly object and its generators
  """
  rp = NCRPoly(pol_gens, bits_exp, field, **kwds)
  return [rp] + rp.gens()

class Poly(dict):
  """
  elements of a multivariate polynomial ring, with bound
  on degrees 0 <= exp_i < self.bits_exp
  """
  def __init__(self,rp,**kw):
    """
    >>> from rmpoly import *
    >>> from fractions import Fraction as mpq
    >>> rp,y,x = rgens('y,x',6,mpq)
    >>> Poly(rp)
    0
    >>> rp(0)
    0
    >>> p1 = x + y**2
    >>> isinstance(p1,Poly)
    True
    """
    self.rp = rp
    dict.__init__(self,**kw)

  def __str__(self):
    """string representation of a polynomial
    Terms are in decreasing order; term coefficients are in front
    >>> from rmpoly import *
    >>> from fractions import Fraction as mpq
    >>> rp,y,x = rgens('y,x',6,mpq)
    >>> x**2/3 + 2*x*y**2/5
     +1/3*x^2 +2/5*x*y^2
    """
    if not self:
      return '0'
    rp = self.rp
    if rp.parens:
      return self.tostr()
    mask_exp = rp.mask_exp
    bits_exp = rp._bits_exp
    pol_gens = rp.pol_gens
    if rp.ordertype == 1:
      mask_vars = rp.mask_vars
    else:
      mask_vars = -1

    s = ''
    a = self.keys()
    a.sort(reverse=True)
    for expv in a:
      c = self[expv]
      if c > 0:
        s += ' +'
      else:
        s += ' -'
      if c < 0:
        c = -c
      if c != 1 or not expv:
        cnt1 = str(c)
      else:
        cnt1 = ''
      # throw away the bits for the hidden variable
      expv = expv & mask_vars
      i = 0
      sa = []
      while(expv):
        exp = expv&mask_exp
        expv >>= bits_exp
        if exp > 1:
          sa.append('%s^%d' % (pol_gens[i],exp))
        if exp == 1:
          sa.append('%s' % pol_gens[i])
        i += 1
      if cnt1:
        sa.append(cnt1)
      sa.reverse()
      s += '*'.join(sa)
    return s

  def tostr(self):
    rp = self.rp
    mask_exp = rp.mask_exp
    bits_exp = rp._bits_exp
    pol_gens = rp.pol_gens
    if rp.ordertype == 1:
      mask_vars = rp.mask_vars
    else:
      mask_vars = -1

    s = ''
    a = self.keys()
    a.sort(reverse=True)
    for expv in a:
      c = self[expv]
      s += ' +(%s)' % c
      if expv:
        s += '*'
      # throw away the bits for the hidden variable
      expv = expv & mask_vars
      i = 0
      sa = []
      while(expv):
        exp = expv&mask_exp
        expv >>= bits_exp
        if exp > 1:
          sa.append('%s^%d' % (pol_gens[i],exp))
        if exp == 1:
          sa.append('%s' % pol_gens[i])
        i += 1
      sa.reverse()
      s += '*'.join(sa)
    return s

  __repr__ = __str__

  def strip(p):
    """eliminate monomials with zero coefficient
    """
    z = p.rp.field(0)
    for k,v in p.items():
      if v == z:
        del p[k]

  def variables(self):
    """return the tuple of the indices of the variables occurring
    in the polynomial p
    >>> from rmpoly import *
    >>> from fractions import Fraction as mpq
    >>> rp,y,x = rgens('y,x',6,mpq)
    >>> p = x**2/3 + 2*x*y**2/5
    >>> p.variables()
    [0, 1]
    >>> rp = RPoly(['x%d' % i for i in range(30)],6,mpq)
    >>> p = rp('x0*x1 + x1*x4 + x4*x8')          
    >>> p.variables()
    [0, 1, 4, 8]
    """
    mask_exp = self.rp.mask_exp
    bits_exp = self.rp._bits_exp
    if self.rp.ordertype:
      mask_vars = self.rp.mask_vars
    else:
      mask_vars = -1
    a = []
    for expv in self:
      if(self[expv] == 0):
        continue
      i = 0
      expv = expv & mask_vars
      while(expv):
        exp = expv&mask_exp
        expv >>= bits_exp
        if exp and i not in a:
          a.append(i)
        i += 1
    a.sort()
    return a

  def __eq__(p1,p2):
    if not p2 or p2 == 0:
      return not p1
    if isinstance(p2,Poly):
      if p1.rp == p2.rp:
        return dict.__eq__(p1,p2)
    else:
      # assume p2 is a coefficient
      if 0 not in p1 or len(p1) > 1:
        return False
      return p1[0] == p1.rp.field(p2)

  def __add__(p1, p2):
    """add two polynomials
    >>> from rmpoly import *
    >>> from fractions import Fraction as mpq
    >>> rp,y,x = rgens('y,x',6,mpq)
    >>> p1 = x + y**2
    >>> p2 = x**2
    >>> p1+p2
     +x^2 +x +y^2
    """
    if not p2:
      return p1.copy()
    if isinstance(p2,Poly):
      if p1.rp == p2.rp:
        p = Poly(p1.rp) 
        for k,v in p1.iteritems():
          if k in p2:
            r = v + p2[k]
            if r:
              p[k] = r
          else:
            p[k] = v
        for k,v in p2.iteritems():
          if k not in p1:
            p[k] = v
        return p
      elif p2.rp == p1.rp.field:
        p = p1.copy()
        if 0 not in p1.keys():
          p[0] = p1.rp.field(p2)
        else:
          if p2 == -p[0]:
            del p[0]
          else:
            p[0] += p2
        return p
      elif p1.rp == p2.rp.field:
        return p2+p1
      else:
        raise ValueError('cannot sum p1 and p2')
    # assume p2 in a number
    else:
      p = p1.copy()
      cp2 = p1.rp.field(p2)
      if not cp2:
        return p
      if 0 not in p1.keys():
        p[0] = cp2
      else:
        if p2 == -p[0]:
          del p[0]
        else:
          p[0] += cp2
      return p

  def __radd__(p1,n):
    # assume n is in p1.rp.field
    p = p1.copy()
    if not n:
      return p
    if 0 not in p1.keys():
      p[0] = p1.rp.field(n)
    else:
      if n == -p[0]:
        del p[0]
      else:
        p[0] += n
    return p

  def __neg__(self):
    return self*(-1)

  def __pos__(self):
    return self

  def __iadd__(p1,p2):
    """inplace add polynomials
    >>> from rmpoly import *
    >>> from fractions import Fraction as mpq
    >>> rp,y,x = rgens('y,x',6,mpq)
    >>> p1 = x + y**2
    >>> p2 = x**2
    >>> p1 += p2
    >>> p1
     +x^2 +x +y^2
    """
    if not p2:
      return p1
    if isinstance(p2,Poly):
      if p1.rp != p2.rp:
        raise ValueError('p1 and p2 must have the same rp')
      dl = []
      for k,v in p1.iteritems():
        if k in p2:
          if p2[k] == -v:
            dl.append(k)
          else:
            p1[k] = v + p2[k]
      for k in p2:
        if k not in p1:
          p1[k] = p2[k]
      for k in dl:
        del p1[k]
      return p1
    else:
      if 0 not in p1.keys():
        p1[0] = p1.rp.field(p2)
      else:
        if p2 == -p1[0]:
          del p1[0]
        else:
          p1[0] += p2
      return p1

  def __sub__(p1,p2):
    if not p2:
      return p1.copy()
    p = Poly(p1.rp)
    if isinstance(p2,Poly):
      if p1.rp == p2.rp:
        for k in p1:
          if k in p2:
            r = p1[k] - p2[k]
            if r:
              p[k] = r
          else:
            p[k] = p1[k]
        for k in p2:
          if k not in p1:
            p[k] = -p2[k]
        return p
      elif p2.rp == p1.rp.field:
        p = p1.copy()
        if 0 not in p1.keys():
           p[0] = -p1.rp.field(p2)
        else:
          if p2 == p[0]:
            del p[0]
          else:
            p[0] -= p2
        return p
      else:
        raise ValueError('cannot coerce p2')
    # assume p2 in a number
    else:
      p2 = p1.rp.field(p2)
      p = copy(p1)
      if 0 not in p1.keys():
        p[0] = -p2
      else:
        if p2 == p[0]:
          del p[0]
        else:
          p[0] -= p2
      return p

  def __rsub__(p1,n):
    # assume n in a number
    p = Poly(p1.rp)
    for expv in p1:
      p[expv] = -p1[expv]
    p += n
    return p

  def __isub__(p1,p2):
    if isinstance(p2,Poly):
      if p1.rp != p2.rp:
        raise ValueError('p1 and p2 must have the same rp')
      dl = []
      for k in p1:
        if k in p2:
          if p2[k] == p1[k]:
            dl.append(k)
          else:
            p1[k] = p1[k] - p2[k]
      for k in p2:
        if k not in p1:
          p1[k] = -p2[k]
      for k in dl:
        del p1[k]
      return p1
    else:
      if 0 not in p1:
        p1[0] = -p2
      else:
        if p1[0] == p2:
          del p1[0]
        else:
          p1[0] -= p2
      return p1

  def univariate_to_dense(p,h):
    """truncate a one-variable series with no parameters and
    convert it to a dense polynomial representation
    """
    if p.ngens > 1:
      raise ValueError, 'only for one variable case'
    keys = p.keys()
    a = zip(keys,p.values())
    a.sort()
    c = [0]*h
    for x in a:
      c[x[0]] = x[1]
    return c

  def __mul__(p1,p2):
    """multiply two polynomials
    >>> from rmpoly import *
    >>> from fractions import Fraction as mpq
    >>> rp,y,x = rgens('y,x',6,mpq)
    >>> p1 = x + y**2
    >>> p2 = x**2
    >>> p1*p2
     +x^3 +x^2*y^2
    """
    p = Poly(p1.rp)
    if not p2:
      return p
    if isinstance(p2,Poly):
      if p1.rp == p2.rp:
        get = p.get
        p2it = p2.items()
        for exp1,v1 in p1.iteritems():
          for exp2,v2 in p2it:
            exp = exp1 + exp2
            p[exp] = get(exp,0) + v1*v2
        p.strip()
        return p
      if p2.rp != p1.rp.field:
        if p1.rp == p2.rp.field:
          p = Poly(p2.rp)
          for exp2,v2 in p2.iteritems():
            p[exp2] = p1*v2
          return p
        else:
          raise ValueError('p1 and p2 must have the same rp')
    # assume p2 in a number
    for exp1,v1 in p1.iteritems():
      v = v1*p2
      if v:
        p[exp1] = v
    return p

  def __rmul__(p1,p2):
    p = Poly(p1.rp)
    if not isinstance(p2,Poly):
      if not p2:
        return p
      for exp1,v1 in p1.iteritems():
        p[exp1] = p2*v1
    else:
      raise NotImplementedError
    return p

  def mul_iadd(p,p1,p2):
    """p += p1*p2
    >>> from rmpoly import *
    >>> from fractions import Fraction as mpq
    >>> rp,y,x = rgens('y,x',6,mpq)
    >>> p1 = x**2 + x
    >>> p2 = x*y + 1
    >>> p3 = y + 1
    >>> p1.mul_iadd(p2,p3)
    >>> p1
     +x^2 +x*y^2 +x*y +x +y +1
    """
    if isinstance(p1,Poly) and isinstance(p2,Poly):
      if p1.rp != p2.rp:
        raise ValueError('p1 and p2 must have the same rp')
      get = p.get
      for exp1,v1 in p1.iteritems():
        for exp2,v2 in p2.iteritems():
          exp = exp1 + exp2
          p[exp] = get(exp,0) + v1*v2
      p.strip()
    else:
      raise NotImplementedError

  def imul_num(p,c):
    """multiply inplace by a number"""
    if not c:
      p.clear()
      return
    for exp in p:
      p[exp] *= c

  def mul_num(p,c):
    if not c:
      return p.rp('0')
    p1 = p.copy()
    for exp in p:
      p1[exp] = c*p[exp]
    return p1

  def __div__(p1,p2):
    """division by a number or exact division by a polynomial
    >>> from rmpoly import *
    >>> from fractions import Fraction as mpq
    >>> rp,y,x = rgens('y,x',6,mpq)
    >>> p = x + 2*y**2
    >>> p/2
     +1/2*x +y^2
    >>> p = (x+y)*(1+x+y)
    >>> p/(x+y)
     +x +y +1
    """
    if isinstance(p2,Poly):
      if p1.rp == p2.rp:
        q,r = p1.division([p2])
        if r:
          raise ValueError, '__div__ performs only division without remainder'
        return q[0]
      elif p2.rp == p1.rp.field:
        p = Poly(p1.rp)
        # if p is not a constant, not implemented
        if p2.keys() != [0]:
          raise NotImplementedError
        else:
          p2 = p2[0]
      else:
        raise NotImplementedError, 'cannot divide p1 by p2'
    # assume p2 in a number
    p = Poly(p1.rp)
    if not p2:
      raise ZeroDivisionError
    try:
      for exp1,v in p1.iteritems():
        p[exp1] = v/p2
    except:
      if not p1.rp.base_field:
        raise NotImplementedError, 'cannot divide p1 by p2'
      else:
        c = p1.rp.base_field(1)/p2
        for exp1,v in p1.iteritems():
          p[exp1] = v*c
    return p

  def mul_mon(p,a, encoded=0):
    """mulltiply p by coeff*x_i^j
       a = (i,j,coeff)
    >>> from rmpoly import *
    >>> from fractions import Fraction as mpq
    >>> rp,y,x = rgens('y,x',6,mpq)
    >>> p1 = x + y**2
    >>> p1.mul_mon((0,2,3))
     +3*x*y^2 +3*y^4
    """
    rp = p.rp
    order = rp.order
    p1 = Poly(rp)
    bits_exp = rp._bits_exp
    j = a[1]
    if encoded:
      expv1 = a[0]
    else:
      i = a[0]
      nb = i*bits_exp
      expv1 = j<<nb
      if rp.ordertype:
        expv1 += (j<<(rp.ngens*bits_exp))*rp.ordertuple[i]
    c = a[-1]
    for expv, v in p.iteritems():
      p1[expv+expv1] = v*c
    return p1

  def mul_trunc(p1,p2,i,h):
    """truncation of p1*p2
    p1 and p2 polynomials
    If h is an integer, 
    let i is the name of the variable x, or its index;
    neglect in p1*p2 the monomials containing x^k, k >= h
    If h is a tuple or list of integers,
    i_j is the name of the variable x_j, or its index, for each
    i_j in i; h_j is the corresponding power in h
    neglect in p1*p2 the monomials containing x_j^k, k >= i_j
    >>> from rmpoly import *
    >>> from fractions import Fraction as mpq
    >>> rp,y,x = rgens('y,x',6,mpq)
    >>> p1 = (x + y)**2
    >>> p1.mul_trunc(x**2,'x',3)
     +x^2*y^2
    >>> rp,z,y,x = rgens('z,y,x',6,mpq)
    >>> p = (1+x+y+z)**3
    >>> p.mul_trunc(1+x,['x','y'],[2,1])
     +x*z^3 +6*x*z^2 +9*x*z +4*x +z^3 +3*z^2 +3*z +1
    """
    rp = p1.rp
    if p1.rp != p2.rp:
      raise ValueError('p1 and p2 must have the same rp')
    mask_exp = rp.mask_exp
    bits_exp = rp._bits_exp
    if hasattr(h,'__hex__'):
      h = int(h)
      if isinstance(i,str):
        iv = rp.gens_dict[i]
      else:
        iv = i
      mask = mask_exp<<(iv*bits_exp)
      maxi = h<<(iv*bits_exp)
      p = Poly(p1.rp)
      get = p.get
# in the case in which iv is the index of the greatest 
# variable (the hidden variable _t in the case of ordertype=1)
# sort the items with respect to expv
      # FIXME put some condition on length?
      #if iv == rp.ngens + rp.ordertype - 1 and len(p1)*len(p2) > 100:
      if iv == rp.ngens + rp.ordertype - 1:
        items2 = p2.items()
        items2.sort()
        for exp1,v1 in p1.iteritems():
          for exp2,v2 in items2:
            exp = exp1 + exp2
            if exp < maxi:
              p[exp] = get(exp,0) + v1*v2
            else:
              break
      else:
        p2it = p2.items()
        for exp1,v1 in p1.iteritems():
          for exp2,v2 in p2it:
            exp = exp1 + exp2
            if exp&mask >= maxi:
              continue
            p[exp] = get(exp,0) + v1*v2
      p.strip()
      return p
    elif isinstance(h,tuple) or isinstance(h,list):
      iv = i
      hv = h[:]
      iv1 = []
      for i in iv:
        if isinstance(i,str):
          i = rp.gens_dict[i]
        iv1.append(i)

      # maskv[i] = mask_exp<<(iv1[i]*bits_exp)
      # is the mask for the exponent of the i-th variable
      # bmask mask for the bits not occupied by the allowed exponents,
      # apart for the place in bmask for the i=0 variable, 
      # for which there is maskv[0];
      # this optimization works provided iv1[0] is the
      # smallest element of iv1; otherwise the smaller variables
      # are not filtered by the bmask test
      miniv1 = min(iv1)
      if miniv1 != iv1[0]:
        imin = iv1.index(miniv1)
        iv1[0],iv1[imin] = iv1[imin],iv1[0]
        hv[0],hv[imin] = hv[imin],hv[0]
      niv = len(iv)
      bmask = mask_exp<<(iv1[0]*bits_exp)
      # mask_items has items (msk,mxv) = (mask_exp<<iv1b,hv[i]<<iv1b)
      # the test (exp&msk) >= mxv is faster; this trick can be used also 
      # in mul_trunc
      mask_items = []
      for i in range(1,niv):
        bhv = bitcount(hv[i]-1)
        iv1b = iv1[i]*bits_exp
        bmask += (mask_exp^((1<<bhv)-1))<<iv1b
        # if hv[i] is a power of 2 the bmask test is sufficient
        if hv[i] != 1<<bhv:
          mask_items.append((mask_exp<<iv1b,hv[i]<<iv1b))
      maxiv0 = hv[0]<<(iv1[0]*bits_exp)
      p = Poly(p1.rp)
      get = p.get
      p2items = p2.items()
      for exp1,v1 in p1.iteritems():
        for exp2,v2 in p2items:
          exp = exp1 + exp2
          if exp&bmask >= maxiv0:
            continue
          accepted = 1
          for msk,mxv in mask_items:
            if (exp&msk) >= mxv:
              accepted = 0
              break
          if accepted:
            p[exp] = get(exp,0) + v1*v2
      p.strip()
      return p

  def square_trunc(p1,i,h):
    """truncation of p1*p1
    If h is an integer, 
    let i is the name of the variable x, or its index;
    neglect in p1*p1 the monomials containing x^k, k >= h
    If h is a tuple or list of integers,
    i_j is the name of the variable x_j, or its index, for each
    i_j in i; h_j is the corresponding power in h
    neglect in p1*p2 the monomials containing x_j^k, k >= i_j
    """
    rp = p1.rp
    if not rp.commuting:
      return p1.mul_trunc(p1,i,h)
    mask_exp = rp.mask_exp
    bits_exp = rp._bits_exp
    if hasattr(h,'__hex__'):
      h = int(h)
      if isinstance(i,str):
        iv = rp.gens_dict[i]
      else:
        iv = i
      mask = mask_exp<<(iv*bits_exp)
      maxi = h<<(iv*bits_exp)
      p = Poly(p1.rp)
      get = p.get
      #print 'DB1', len(p1)
# in the case in which iv is the index of the greatest 
# variable (the hidden variable _t in the case of ordertype=1)
# sort the items with respect to expv, which implies 
# sorting with respect to the powere of iv; then if
# exp1+exp2 >= maxi  the corresponding monomial has
# power of iv >= h, so it is truncated; for the same
# exp1, the terms subsequent to exp2 correspond to the same
# power or higher, so they are also truncated
# TODO treat in a similar way the case in which iv is not 
# the greatest variable
      if iv == rp.ngens + rp.ordertype - 1:
        items = p1.items()
        items.sort()
        #print 'DB2',items
        for i in xrange(len(items)):
          exp1,v1 = items[i]
          for j in range(i):
            exp2,v2 = items[j]
            exp = exp1 + exp2
            if exp&mask < maxi:
              p[exp] = get(exp,0) + v1*v2
            else:
              break
      else:
        keys = p1.keys()
        for i in range(len(keys)):
          exp1 = keys[i]
          v1 = p1[exp1]
          for j in range(i):
            exp2 = keys[j]
            exp = exp1 + exp2
            if exp&mask < maxi:
              p[exp] = get(exp,0) + v1*p1[exp2]

      p.imul_num(2)
      for k,v in p1.iteritems():
        if (2*k)&mask < maxi:
          p[2*k] = get(2*k,0) + v**2
      p.strip()
      return p

    elif isinstance(h,tuple) or isinstance(h,list):
      iv = i
      hv = h[:]
      iv1 = []
      for i in iv:
        if isinstance(i,str):
          i = rp.gens_dict[i]
        iv1.append(i)

      # maskv[i] = mask_exp<<(iv1[i]*bits_exp)
      # is the mask for the exponent of the i-th variable
      # bmask mask for the bits not occupied by the allowed exponents,
      # apart for the place in bmask for the i=0 variable, 
      # for which there is maskv[0];
      # this optimization works provided iv1[0] is the
      # smallest element of iv1; otherwise the smaller variables
      # are not filtered by the bmask test
      miniv1 = min(iv1)
      if miniv1 != iv1[0]:
        imin = iv1.index(miniv1)
        iv1[0],iv1[imin] = iv1[imin],iv1[0]
        hv[0],hv[imin] = hv[imin],hv[0]
      niv = len(iv)
      bmask = mask_exp<<(iv1[0]*bits_exp)
      # mask_items has items (msk,mxv) = (mask_exp<<iv1b,hv[i]<<iv1b)
      # the test (exp&msk) >= mxv is faster; this trick can be used also 
      # in mul_trunc
      mask_items = []
      for i in range(1,niv):
        bhv = bitcount(hv[i]-1)
        iv1b = iv1[i]*bits_exp
        bmask += (mask_exp^((1<<bhv)-1))<<iv1b
        # if hv[i] is a power of 2 the bmask test is sufficient
        if hv[i] != 1<<bhv:
          #print 'DB1',hv[i], 1<<bhv
          mask_items.append((mask_exp<<iv1b,hv[i]<<iv1b))
      maxiv0 = hv[0]<<(iv1[0]*bits_exp)
      p = Poly(p1.rp)
      #assert p1.rp == p2.rp
      get = p.get
      keys = p1.keys()
      for i in range(len(keys)):
        exp1 = keys[i]
        v1 = p1[exp1]
        for j in range(i):
          exp2 = keys[j]
          exp = exp1 + exp2
          if exp&bmask >= maxiv0:
            continue
          accepted = 1
          for msk,mxv in mask_items:
            if (exp&msk) >= mxv:
              accepted = 0
              break
          if accepted:
            p[exp] = get(exp,0) + v1*p1[exp2]
      p.imul_num(2)
      for k,v in p1.iteritems():
        exp = 2*k
        if exp&bmask >= maxiv0:
          continue
        accepted = 1
        for msk,mxv in mask_items:
          if (exp&msk) >= mxv:
            accepted = 0
            break
        if accepted:
          p[exp] = get(exp,0) + v**2
    p.strip()
    return p



  def pow_trunc(self,n,i,h):
    """truncation of self**n using mul_trunc
    >>> from fractions import Fraction as mpq
    >>> rp,y,x = rgens('y,x',6,mpq)
    >>> p1 = 1 + x + y
    >>> p1.pow_trunc(3,'y',1)
     +x^3 +3*x^2 +3*x +1
    >>> p1.pow_trunc(3,['x','y'],[2,1])
     +3*x +1
    """
    rp = self.rp
    if not hasattr(n,'__hex__'):
      raise NotImplementedError
    n = int(n)
    if n == 0:
      if self:
        return rp('1')
      else:
        raise ValueError
    if n < 0:
      # for univariate series series_inversion is faster than pow)trunc
      if hasattr(h,'__hex__'):
        p1 = self.pow_trunc(-n,i,h)
        return p1.series_inversion(i,h)
      else:
        p1 = self.series_inversion(i,h)
        return p1.pow_trunc(-n,i,h)
    if n == 1:
      return self.trunc(i,h)
    if n == 2:
      return self.square_trunc(i,h)
    if n == 3:
      p2 = self.square_trunc(i,h)
      return p2.mul_trunc(self,i,h)
    p = rp('1')
    while 1:
      if n&1:
        p = self.mul_trunc(p,i,h)
        n -= 1
        if not n:
          break
      self = self.square_trunc(i,h)
      n = n // 2
    return p

  def trunc(p1,i,h):
    """monomials containing x^k, k >= h neglected
    i is the name of the variable x, or its index
    >>> from fractions import Fraction as mpq
    >>> rp,y,x = rgens('y,x',6,mpq)
    >>> p = (1 + x + y)**3
    >>> p.trunc('x',2)
     +3*x*y^2 +6*x*y +3*x +y^3 +3*y^2 +3*y +1
    >>> p.trunc(['x','y'],[2,2])
     +6*x*y +3*x +3*y +1
    """
    rp = p1.rp
    mask_exp = rp.mask_exp
    bits_exp = rp._bits_exp
    if hasattr(h,'__hex__'):
      h = int(h)
      if isinstance(i,str):
        i = rp.gens_dict[i]
      mask = mask_exp<<(i*bits_exp)
      maxi = h<<(i*bits_exp)
      p = Poly(p1.rp)
      for exp1 in p1:
        if exp1&mask >= maxi:
          continue
        p[exp1] = p1[exp1]
      return p
    elif isinstance(h,tuple) or isinstance(h,list):
      iv = i
      hv = h
      iv1 = []
      for i in iv:
        if isinstance(i,str):
          i = rp.gens_dict[i]
        iv1.append(i)
      maskv = []
      maxiv = []
      niv = len(iv)
      for i in range(niv):
        mask = mask_exp<<(iv1[i]*bits_exp)
        maxi = hv[i]<<(iv1[i]*bits_exp)
        maskv.append(mask)
        maxiv.append(maxi)
      p = Poly(p1.rp)
      rangeniv = range(niv)
      for exp in p1:
        accepted = 1
        for i in rangeniv:
          if (exp&maskv[i]) >= maxiv[i]:
            accepted = 0
            break
        if not accepted:
          continue
        p[exp] = p1[exp]
      return p

  def iadd_mon(self,a,encoded=0):
    """dd the monomial coeff*x0^i0*x1^i1*...
       a = ((i0,i1,...),coeff)
    """
    coeff = a[1]
    if encoded:
      expv = a[0]
    else:
      expv = self.rp.code_expv(a[0])
    if expv in self:
      self[expv] += coeff
      if self[expv] == 0:
        del self[expv]
    else:
      self[expv] = coeff

  def total_degree(self):
    """maximum of the degrees of the terms in the polynomial
    >>> from fractions import Fraction as mpq
    >>> rp,y,x = rgens('y,x',6,mpq)
    >>> p = 1 + x + x*y
    >>> p.total_degree()
    2
    """
    if not self:
      return 0
    rp = self.rp
    if rp.order == 'grlex':
      expvmax = max(self)
      deg = expvmax>>(rp.ngens*rp._bits_exp)
    else:
      deg = 0
      for expv in self:
        deg1 = sum(rp.decode_expv(expv))
        deg = max(deg,deg1)
    return deg

  def derivative_t(self,n):
    """
    return the sum of the monomials of p, each multiplied by its
    total degree if n = 1, else divided by the total degree

    _t*d/d _t p for p with grlex order  for n == 1
    integral d _t (p/_t)                for n != 1 (p without constant terms)
    """
    rp = self.rp
    if rp.order != 'grlex':
      return NotImplementedError
    nb = rp.ngens*rp._bits_exp
    p1 = Poly(rp)
    if n == 1:
      for expv in self:
        if expv:
          exp = expv >> nb
          p1[expv] = self[expv]*exp
    else:
      if 0 in self:
        raise ValueError, 'no constant term allowed'
      try:
        for expv in self:
          exp = expv >> nb
          p1[expv] = self[expv]/exp
      except:
        if rp.base_field:
          one = rp.base_field(1)
          for expv in self:
            exp = expv >> nb
            print 'DB1', exp, one/exp
            p1[expv] = self[expv]*(one/exp)
    return p1

  def derivative(self, n):
    """derivative of p with respect to x_n; the argument n is the
    index of the variable x_n = self.rp.gens()[n], or the 
    corresponding string.
    >>> from fractions import Fraction as mpq
    >>> rp,y,x = rgens('y,x',6,mpq)
    >>> p = x + x**2*y**3
    >>> p.derivative('x')
     +2*x*y^3 +1
    """
    rp = self.rp
    mask_exp = rp.mask_exp
    bits_exp = rp._bits_exp
    pol_gens = rp.pol_gens
    if n == '_t':
      return self.derivative_t()
    if n.__class__ == str:
      n = rp.gens_dict[n]
    else:
      n = int(n)
    if rp.ordertype:
      expv_ddeg = (1<<(rp.ngens*bits_exp))*rp.ordertuple[n]
    else:
      expv_ddeg = 0
    if n.__class__ == str:
      n = rp.gens_dict[n]
    p1 = Poly(rp)
    nb = n*bits_exp
    snb = 1<<nb
    # mask to check presence of variable n 
    mask = mask_exp<<nb
    for expv in self:
      if expv&mask:
        exp = (expv>>nb)&mask_exp
        p1[expv-snb-expv_ddeg] = self[expv]*exp
    return p1


  def integrate(self, n):
    """ integrate p with respect to x_n; the argument n is the
    index of the variable x_n = self.rp.gens()[n], or the 
    corresponding string.
    >>> from fractions import Fraction as mpq
    >>> rp,y,x = rgens('y,x',6,mpq)
    >>> p = x + x**2*y**3
    >>> p.integrate('x')
     +1/3*x^3*y^3 +1/2*x^2
    """
    mask_exp = self.rp.mask_exp
    bits_exp = self.rp._bits_exp
    pol_gens = self.rp.pol_gens
    if n.__class__ == str:
      n = self.rp.gens_dict[n]
    if self.rp.ordertype:
      expv_ddeg = (1<<(self.rp.ngens*bits_exp))*self.rp.ordertuple[n]
    else:
      expv_ddeg = 0
    p1 = Poly(self.rp)
    nb = n*bits_exp
    snb = 1<<nb
    mask = mask_exp<<nb
    for expv in self:
      exp = (expv>>nb)&mask_exp
      p1[expv+snb+expv_ddeg] = self[expv]/(exp+1)
    return p1


  def copy(self):
    return copy(self)

  def coefficient_t(self,t):
    """coefficient x_i^j of p
    t=(i,j) corresponds to x_i^j
    for j=0 it gives the terms independent from the variable x_i
    >>> from fractions import Fraction as mpq
    >>> rp,y,x = rgens('y,x',6,mpq)
    >>> p = (1+x+y)**3
    >>> p.coefficient_t((0,2))
     +3*x +3
    >>> p.coefficient_t((0,0))
     +x^3 +3*x^2 +3*x +1
    """
    rp = self.rp
    mask_exp = rp.mask_exp
    bits_exp = rp._bits_exp
    i,j = t
    nb = i*bits_exp
    mask1 = mask_exp<<nb
    expv1 = j<<nb
    p = Poly(rp)
    if rp.ordertype:
      expv1 += (j<<(rp.ngens*bits_exp))*rp.ordertuple[i]
      # mask_vars used to exclude the exponent of the hidden variable
      mask_vars = rp.mask_vars
      expv1a=expv1&mask_vars
      for expv in self:
        if (expv&mask1)&mask_vars == expv1a:
          p[expv-expv1] = self[expv]
    else:
      for expv in self:
        if expv&mask1 == expv1:
          p[expv-expv1] = self[expv]
    return p

  def coefficient(self,p1):
    """the coefficient of a monomial p1 is the sum ot the terms in 
    self which have the same degrees in the variables present in p1,
    divided by p1
    >>> from fractions import Fraction as mpq
    >>> rp,y,x = rgens('y,x',6,mpq)
    >>> p = (1 + x + y)**3
    >>> p.coefficient(y**2)
     +3*x +3
    """
    rp = self.rp
    mask_exp = rp.mask_exp
    bits_exp = rp._bits_exp
    order = rp.order
    k = p1.keys()
    if len(k) != 1:
      raise TypeError('the argument of coeff must be a monomial')
    expv1 = k[0]
    # mask1 used to select the exponents of the variables present in p1
    v1 = p1.variables()
    mask1 = 0
    for i in v1:
      mask1 += mask_exp<<(i*bits_exp)
    p = Poly(rp)
    if rp.ordertype:
      # mask_vars used to exclude the exponent of the hidden variable
      mask_vars = rp.mask_vars
      expv1a=expv1&mask_vars
      for expv in self:
        if (expv&mask1)&mask_vars == expv1a:
          p[expv-expv1] = self[expv]
    else:
      for expv in self:
        if expv&mask1 == expv1:
          p[expv-expv1] = self[expv]
    return p

  def coeff(self,p1):
    """monomial coefficient
    >>> from fractions import Fraction as mpq
    >>> rp,y,x = rgens('y,x',6,mpq)
    >>> p = (1 + x + y)**3
    >>> print p.coeff(y**2)
    3
    >>> print p.coeff(x*y)
    6
    """
    expv = p1.leading_expv()
    if expv in self:
      return self[expv]
    else:
      return 0


  def expv_mask(self):
    """given a monomial, return expv1,mask1
    expv1 encoded exponents
    mask1 mask for the bits of the exponents of the
    variables present in the monomial
    """
    rp = self.rp
    mask_exp = rp.mask_exp
    bits_exp = rp.bits_exp
    k = self.keys()
    if len(k) != 1:
      raise TypeError('the argument of coeff must be a monomial')
    expv1 = k[0]
    v1 = self.variables()
    mask1 = 0
    for i in v1:
      mask1 += mask_exp<<(i*bits_exp)
    return expv1, mask1

  def coefficient1(self,expv1,mask1):
    """coefficient of a monomial, fast version;
    expv1 encoded exponents
    mask1 mask for the bits of the exponents of the
    variables present in the monomial
    """
    rp = self.rp
    mask_exp = rp.mask_exp
    bits_exp = rp.bits_exp
    p = Poly(rp)
    for expv in self:
      if expv&mask1 == expv1:
        p[expv-expv1] = self[expv]
    return p

  def iadd_m_mul_q(p1,p2,(expva,c)):
    """ p1 += p2*m
     m = monomial with m[expva] = c
    """
    get = p1.get
    for k,v in p2.iteritems():
      ka = k+expva
      p1[ka] = get(ka,0) + v*c
    p1.strip()
    return p1

  def leading_expv(self):
    """leading expv 
    """
    if self:
      #return max(self.keys())
      return max(self)
    else:
      return None

  def leading_term(self):
    """leading term according to grlex ordering with reversed order
    of variables
    >>> from fractions import Fraction as mpq
    >>> rp,y,x = rgens('y,x',6,mpq)
    >>> p = x**2 + y**3
    >>> p.leading_term()
     +x^2
    """
    p = Poly(self.rp)
    expv = self.leading_expv()
    if expv == None:
      return self.rp('0')
    p[expv] = self[expv]
    return p

  def division(self, fv):
    """division algorithm, see [CLO] p64
    fv array of polynomials
       return qv, r such that
       self = sum(fv[i]*qv[i]) + r
    >>> from fractions import Fraction as mpq
    >>> rp,y,x = rgens('y,x',6,mpq)
    >>> p = (1+x)*(1+y) + 2*(1+x)*(1+y)**2
    >>> p.division([1+x,1+y])
    ([ +2*y^2 +5*y +3, 0], 0)
    >>> p.division([x,y])
    ([ +2*y^2 +5*y +3,  +2*y +5],  +3)
    """
    rp = self.rp
    if not self:
      return Poly(rp), Poly(rp)
    for f in fv:
      if f.rp != rp:
        raise ValueError('self and f must have the same rp')
    s = len(fv)
    qv = [Poly(rp) for i in range(s)]
    p = self.copy()
    r = Poly(rp)
    maskof = rp.maskof
    order = self.rp.order
    expvs = [fx.leading_expv() for fx in fv]
    while p:
      i = 0
      divoccurred = 0
      while i < s and divoccurred == 0:
        expv = p.leading_expv()
        expv1 = expv - expvs[i]
        # test if the leading term of fv[i] divides the leading term of p 
        if expv1 >= 0 and expv1&maskof == 0:
          c = p[expv]/fv[i][expvs[i]]
          qv[i].iadd_mon((expv1,c),1)
          p.iadd_m_mul_q(fv[i],(expv1,-c))
          # in grlex order the powers of the variables is bounded
          # by the total degree of self; in other orders this is 
          # not generally true;
          # check that there is no overflow in these cases
          if order != 'grlex':
            bov = 0
            for expv2 in p:
              if expv2&maskof:
                bov = 1
                break
            if bov:
              raise RPolyOverflowError('Division')

          divoccurred = 1
          #if expv == 0:
          #  print 'DB5 expv==0'
          #  break
        else:
          i += 1
      if not divoccurred:
        expv =  p.leading_expv()
        #print 'DB3a p=',p,expv, p.leading_expv()
        r.iadd_mon((expv,p[expv]),1)
        del p[expv]
    if expv == 0:
      r += p
    return qv, r

  def mod1(self, fv):
    """fv array of (expv, p), p monic
       return r such that
       self = sum(fv[i]*qv[i]) + r
    """
    rp = self.rp
    if not self:
      return Poly(rp)
    s = len(fv)
    p = self.copy()
    r = Poly(rp)
    maskof = rp.maskof
    order = self.rp.order
    expvs = [fx[0] for fx in fv]
    fv = [fx[1] for fx in fv]
    expv = p.leading_expv()
    while p:
      i = 0
      divoccurred = 0
      while i < s and divoccurred == 0:
        expv1 = expv - expvs[i]
        # test if the leading term of fv[i] divides the leading term of p 
        if expv1 >= 0 and expv1&maskof == 0:
          p.iadd_m_mul_q(fv[i],(expv1,-p[expv]))
          expv = p.leading_expv()
          if order != 'grlex':
            bov = 0
            for expv2 in p:
              if expv2&maskof:
                bov = 1
                break
            if bov:
              raise RPolyOverflowError('Division with lex')

          divoccurred = 1
        else:
          i += 1
      if not divoccurred:
        r.iadd_mon((expv,p[expv]),1)
        del p[expv]
        expv = p.leading_expv()
    if expv == 0:
      r += p
    return r

  def reduce(self, fv):
    """rest in the division algorithm
    >>> from fractions import Fraction as mpq
    >>> rp,y,x = rgens('y,x',6,mpq)
    >>> p = (1+x)*(1+y) + 2*(1+x)*(1+y)**2
    >>> p.reduce([x,y])
     +3
    """
    return self.division(fv)[1]

  def square(p1):
    rp = p1.rp
    if not rp.commuting:
      return p1*p1
    p = Poly(rp)
    get = p.get
    keys = p1.keys()
    for i in range(len(keys)):
      k1 = keys[i]
      pk = p1[k1]
      for j in range(i):
        k2 = keys[j]
        exp = k1+k2
        p[exp] = get(exp,0) + pk*p1[k2]
    p.imul_num(2)
    for k,v in p1.iteritems():
      p[2*k] = get(2*k,0) + v**2
    p.strip()
    return p

  def miller_pow(self, n, mc=None):
    """power of a polynomial using Miller recursive algorithm
    n  integer power
    mc multinomial coefficients; typically it takes longer to compute
    them than to compute the polynomial, so it can be convenient to cache 
    them.
    Adapted from sympycore-0.1 by Pearu Peterson
    Copyright (c) 2007-2008 Pearu Peterson
    """
    rp = self.rp
    if not rp.commuting:
      return p1**n
    if n == 0:
      return rp('1')
    if n == 1:
      return copy(self)
    if hasattr(n,'__hex__'):
      n = int(n)
    else:
      raise NotImplementedError
    if len(self) == 1:
      p = Poly(rp)
      k, v = self.items()[0]
      # treat case abs(v) = 1 separately to deal with the case in which
      # n is too large to be allowed in v**n
      if v == 1:
        p[k*n] = v
      elif v == -1:
        if n%2 == 0:
          p[k*n] = -v
        else:
          p[k*n] = v
      else:
        p[k*n] = v**n
      return p
    #if n == 2:
    #  return self.square()
    items = self.items()
    p = Poly(self.rp)
    get = p.get
    if not mc:
      m = len(self)
      mc = multinomial_coefficients(m,n)
    for k,c_kn in mc.iteritems():
      exp = 0
      v = c_kn
      for i,e in enumerate(k):
        if e:
          exps, coeff = items[i]
          exp += exps*e
          v *= coeff ** e
      p[exp] = get(exp,0) + v
    p.strip()
    return p

  def __pow__(self,n):
    """power of a polynomial
    >>> from fractions import Fraction as mpq
    >>> rp,y,x = rgens('y,x',6,mpq)
    >>> assert (x+y)**4 == rp('x^4 +4*x^3*y +6*x^2*y^2 +4*x*y^3 +y^4')
    """
    rp = self.rp
    # test if n is an integer
    if not hasattr(n,'__hex__'):
      raise NotImplementedError
    n = int(n)
    if n < 0:
      raise ValueError('n >= 0 is required')
    if n == 0:
      if self:
        return rp('1')
      else:
        raise ValueError
    elif len(self) == 1:
      p = Poly(rp)
      k, v = self.items()[0]
      # treat case abs(v) = 1 separately to deal with the case in which
      # n is too large to be allowed in v**n
      if v == 1:
        p[k*n] = v
      elif v == -1:
        if n%2 == 0:
          p[k*n] = -v
        else:
          p[k*n] = v
      else:
        p[k*n] = v**n
      return p
    elif n == 1:
      return copy(self)
    elif n == 2:
      return self.square()
    elif n == 3:
      return self*self.square()
    p = rp('1')
    while 1:
      if n&1:
        p = p*self
        n -= 1
        if not n:
          break
      self = self.square()
      n = n // 2
    return p

  def has_constant_term(p,iv):
    rp = p.rp
    if isinstance(iv, basestring):
      if iv == '_t':
        if rp.order == 'grlex':
          return 0 in p
        else:
          # TODO case of other weight ordering
          raise NotImplementedError
      iv = [rp.gens_dict[iv]]
    else:
      if isinstance(iv[0], basestring):
        d = rp.gens_dict
        iv = [d[x] for x in iv]
    mask_exp = rp.mask_exp
    bits_exp = rp._bits_exp
    mask = 0
    for i in iv:
      mask += mask_exp<<(i*bits_exp)
    for expv in p:
      if not mask&expv:
        return True
    return False

  def _series_inversion1(p,iv,nv):
    """univariate series inversion 1/p
    iv name of the series variable
    nv precision of the series

    The Newton method is used.
    """
    rp = p.rp
    if 0 not in p:
      raise ValueError('no constant term in series')
    if (p - p[0]).has_constant_term(iv):
      raise ValueError, 'p cannot contain a constant term depending on parameters'
    if p[0] != rp.field(1):
      p1 = rp(1)/p[0]
    else:
      p1 = rp(1)
    #print 'DB', p, iv, nv
    for prec in newton_method_sizes(nv):
        #print 'DB2', prec
        tmp = p1.square()
        tmp = tmp.mul_trunc(p, iv, prec)
        #print 'DB3', prec, p1-tmp 
        p1 = 2*p1 - tmp
    return p1

  def _series_inversion_nc(p,iv,nv):
    """
    a = A^-1, c =  B*a
    (A + B)^-1 = a*(1 - c + c^2 - ...)
    """
    rp = p.rp
    if p[0] != rp.field(1):
      a = p[0]**-1
    else:
      a = rp.field(1)
    b = p - p[0]
    c = -b*a
    cn = rp(1)
    r = rp(1)
    z = rp(0)
    while 1:
      cn = cn.mul_trunc(c, iv,nv)
      if cn == z:
        break
      r += cn
    #r = r.trunc(iv,nv)
    r = a*r
    return r



  def series_inversion(p,iv,nv):
    """multivariate series inversion 1/p
    iv list of variable names or variable indices
    nv list of truncations for these variables
    In the case of one variable it can also be:
      iv variable name or variable index (0)
      nv truncation integer for the variable
    p is a series with O(x_1^n_1*..x_m^n_m) in 
    variables x_k with index or name iv[k-1]
    p has constant term different from zero
    >>> from rmpoly import *
    >>> from fractions import Fraction as mpq
    >>> rp,a,x = rgens('a,x',6,mpq)
    >>> (1 + a*x).series_inversion('x', 6)
     -x^5*a^5 +x^4*a^4 -x^3*a^3 +x^2*a^2 -x*a +1
    """
    rp = p.rp
    if 0 not in p:
      raise ValueError('no constant term in series')
    if (p-p[0]).has_constant_term(iv):
      raise NotImplementedError, 'p-p[0] must not have a constant term in the series variables'
    if not rp.commuting:
      return p._series_inversion_nc(iv,nv)
    if hasattr(nv,'__hex__'):
      return p._series_inversion1(iv, nv)
    iv1 = []
    for i in iv:
      if isinstance(i,str):
        i = rp.gens_dict[i]
      iv1.append(i)
    if 0 not in p:
      raise ValueError('no constant term in series')
    a0 = p[0]
    a = 1 - p/a0
    # TODO benchmark
    if len(p) < 100:
      an = rp('1')
      r = rp('1')
      while 1:
        an = an.mul_trunc(a, iv1,nv)
        if not an:
          break
        r += an
      #r = r.trunc(iv,nv)
      r = r/a0
      return r
    else:
      # tn = A + A^2 + ... + A^n
      # t2n = tn + A^n*tn
      t = a.copy()
      p = a
      #i = 1
      while 1:
        t += t.mul_trunc(p,iv1,nv)
        p = p.square_trunc(iv1,nv)
        if not p:
          break
      t += 1
      t = t/a0
      return t

  def series_reversion(p,i,n,j):
    """reversion of a series
    p is a series with O(x^n) of the form 
    p = a*x + f(x) 
    where a is a number different from 0
    f(x) = sum( a_k*x_k, k in range(2,n))
    a_k can depend polynomially on other variables, not indicated.
    x variable with index i, or with name i
    y variable with index j, or with name j
    solve p = y, that is given
    a*x + f(x) - y = 0
    find the solution x = r(y) up to O(y^n)
  
    if r_i is the solution at order i
    a*r_i + f(r_i) - y = O(y^(i+1))
    r_(i+1) = r_i + e such that
    a*r_(i+1) + f(r_(i+1)) - y = O(y^(i+2))
    a*e + f(r_i) = O(y^(i+2))
    e = -f(r_i)/a
    so that one uses the recursion relation
    r_(i+1) = r_i -f(r_i)/a
    with the boundary condition
    r_1 = y
    >>> from rmpoly import *
    >>> from fractions import Fraction as mpq
    >>> rp,y,x = rgens('y,x',6,mpq)
    >>> p = x+x**2
    >>> p1 = p.series_reversion('x',4,'y')
    >>> p1
     +2*y^3 -y^2 +y
    >>> (p1 + p1**2).trunc('y',4)
     +y
    """
    rp = p.rp
    if not rp.commuting:
      raise NotImplementedError
    if i.__class__ == str:
      nx = i
      i = p.rp.gens_dict[i]
    else:
      nx = rp.pol_gens[i]
    y = rp(j)
    if j.__class__ == str:
      ny = j
      j = p.rp.gens_dict[j]
    else:
      ny = rp.pol_gens[j]
    if p.coefficient_t((i,0)):
      raise ValueError('part independent from %s must be 0' % nx)
    a = p.coefficient_t((i,1))
    assert 0 in a and len(a) == 1
    a = a[0]
    r = rp(ny)/a
    for i in range(2,n):
      sb = Subs(rp,rp,{nx:r})
      sp = sb.subs_trunc(p,ny,i+1)
      sp = sp.coefficient_t((j,i))*y**i
      r -= sp/a
    return r

  def series_reversion2(p,i1,n1,iv,nv,j):
    """reversion of a multivariate series
    p is a multivariate series with O(x^n1*x_1^n_1*..*x_m^n_m)
    of the form 
    p = a(x_1,..,x_n)*x + f(x,x_1,..,x_n) 
    with a(0,..,0) number different from 0
    and 
    f(x,x_1,..,x_m) = sum( a_k*x1_k, k in range(2,n1))
    a, a_k  series in x_1,..,x_m with O(x_1^n_1*..*x_m^n_m) 
    a_k can depend polynomially on other variables, not indicated.
    x variable with index i1, or with name i1; n=n1
    x_k variable with index or name iv[k-1]
    (In the case m=1 ine can integers iv, nv instead of tuples)
    y variable with index j, or with name j
    solve p = y, that is given
    a(x_1,.,x_m)*x + f(x,x_1,..,x_m) - y = 0
    find the solution x = r(x_1,..,x_m,y) up to O(y^n1*x_1^n_1*...x_m^n_m)
  
    a1 = 1/a series inversion
    if r_i is the solution at order i
    a*r_i + f(r_i) - y = O(y^(i+1))
    r_(i+1) = r_i + e such that
    a*r_(i+1) + f(r_(i+1)) - y = O(y^(i+2)*x_1^n_1*...x_m^n_m)
    a*e + f(r_i) = O(y^(i+2)*x_1^n_1*...x_m^n_m)
    e = -f(r_i)*a1
    r_1 = y
    r_(i+1) = r_i -f(r_i)*a1
    >>> from rmpoly import *
    >>> from fractions import Fraction as mpq
    >>> rp,z,y,x = rgens('z,y,x',6,mpq)
    >>> p = x*(1+y) + x**2*(1+y)
    >>> p.series_reversion2('x',2,'y',3,'z')
     +y^2*z -y*z +z
    >>> p1 = p.series_reversion2('x',2,'y',3,'z')
    >>> p1
     +y^2*z -y*z +z
    >>> (p1*(1+y) + p1**2*(1+y)).trunc('z',2).trunc('y',3)
     +z
    >>> rp,w,z,y,x = rgens('w,z,y,x',6,mpq)
    >>> p = x*(1+y+z) + x**2*(1+y+z**2)
    >>> p.series_reversion2('x',2,['y','z'],[3,3],'w')
     +6*y^2*z^2*w -3*y^2*z*w +y^2*w -3*y*z^2*w +2*y*z*w -y*w +z^2*w -z*w +w
    >>> p1 = p.series_reversion2('x',2,['y','z'],[3,3],'w')
    >>> p1
     +6*y^2*z^2*w -3*y^2*z*w +y^2*w -3*y*z^2*w +2*y*z*w -y*w +z^2*w -z*w +w
    >>> Subs(rp,rp,{'x':p1}).subs_trunc(p,['w','y','z'],[2,3,3])
     +w
    """
    rp = p.rp
    if not rp.commuting:
      raise NotImplementedError
    if i1.__class__ == str:
      nx1 = i1
      i1 = p.rp.gens_dict[i1]
    else:
      nx1 = rp.pol_gens[i1]
    # test if nv is an integer
    if hasattr(nv,'__hex__'):
      iv = [iv]
      nv = [nv]
    n2 = len(nv)
    nv = map(int,nv)
    iv1 = []
    for i in iv:
      if isinstance(i,str):
        i = rp.gens_dict[i]
      iv1.append(i)
    y = rp(j)
    if j.__class__ == str:
      ny = j
      j = p.rp.gens_dict[j]
    else:
      ny = rp.pol_gens[j]
    if p.coefficient_t((i1,0)):
      raise ValueError('part independent from %s must be 0' % nx1)
    a = p.coefficient_t((i1,1))
    a1 = a.series_inversion(iv1,nv)
    r = rp(ny)*a1
    for i in range(2,n1):
      sb = Subs(rp,rp,{nx1:r})
      sp = sb.subs_trunc(p,[ny]+iv1,[i+1]+nv)
      sp = sp.coefficient_t((j,i))*y**i
      r -= sp.mul_trunc(a1,[ny]+iv1,[i+1]+nv)
    return r

  def series_from_list(p, c, iv,nv,concur=1):
    """series sum c[n]*p^n
    reduce the number of multiplication summing concurrently 
    ax = [1,p,p^2,..,p^(J-1)]
    s = sum(c[i]*ax[i] for i in range(0,J)) +
        sum(c[i]*ax[i] for i in range(J,2*J))*p^J +
        sum(c[i]*ax[i] for i in range(2*J,3*J))*p^(2*J) + ...+
        sum(c[i]*ax[i] for i in range((K-1)*J,K*J))*p^((K-1)*J)
    with K >= (n+1)/J
    """
    rp = p.rp
    #print 'DB10', dict(p)
    assert 0 not in p
    n = len(c)
    if not concur:
      q = rp('1')
      s = c[0]*q
      for i in range(1,n):
        q = q.mul_trunc(p,iv,nv)
        s += c[i]*q
      return s
    J = int(math.sqrt(n) + 1)
    if n%J == 0:
      K = n/J
    else:
      K = int(n/J) + 1
    ax = [rp('1')]
    b = 1
    q = rp('1')
    # TODO get a good value
    if len(p) < 20:
      for i in range(1,J):
        q = q.mul_trunc(p, iv,nv)
        ax.append(q)
    else:
      for i in range(1,J):
        if i%2 == 0:
          q = ax[i//2].square_trunc(iv,nv)
        else:
          q = q.mul_trunc(p, iv,nv)
        ax.append(q)
    # optimize using square_trunc
    pj = ax[-1].mul_trunc(p, iv,nv)
    b = rp('1')
    s = rp('0')
    for k in range(K-1):
      r = J*k
      s1 = c[r]
      for j in range(1,J):
        s1 += c[r+j]*ax[j]
      s1 = s1.mul_trunc(b,iv,nv)
      s += s1
      b = b.mul_trunc(pj,iv,nv)
      if not b:
        break
    k = K-1
    r = J*k
    if r < n:
      s1 = c[r]*rp('1')
      for j in range(1,J):
        if r+j >= n:
          break
        s1 += c[r+j]*ax[j]
      s1 = s1.mul_trunc(b,iv,nv)
      s += s1
    return s

  def fun(p,f,*args):
    """
    function of a multivariate series computed by substitution

    p multivariate series
    f method name or function
    args[:-2] arguments of f, apart from the first one
    args[-2] = iv names of the series variables
    args[-1] = nv list of the precisions of the series variables
    The case with f method name is used to compute tan and nth_root
    of a multivariate series:
    p.fun('tan',iv,nv)
    compute _x.tan(iv,sum(nv)), then substitute _x with p
    p.fun('nth_root',n,iv,nv)
    compute (_x + p[0]).nth_root(n,'_x',sum(nv)), then substitute _x
    with p - p[0]
    example with f function:
    f = Poly.exp_series0
    p1 = p.fun(Poly.exp_series0,['y','x'],[h,h])
    """
    rp = p.rp
    rp1 = RPoly(['_x'], rp.bits_exp,rp.field)
    _x = rp1.gens()[0]
    h = args[-1]
    if not hasattr(h,'__hex__'):
      h = sum(args[-1])
    args1 = args[:-2] + ('_x', h)
    # separate the constant term of the series
    # compute the univariate series f(_x,..,'x',sum(nv))
    # or _x.f(...,'x',sum(nv)
    if 0 in p:
      x1 = _x + p[0]
      p1 = p - p[0]
    else:
      x1 = _x
      p1 = p
    if isinstance(f,basestring):
      q = getattr(x1,f)(*args1)
    else:
      q = f(x1,*args1)
    a = zip(q.keys(),q.values())
    a.sort()
    c = [0]*h
    for x in a:
      c[x[0]] = x[1]
    p1 = p1.series_from_list(c,args[-2],args[-1])
    return p1


  def subs(p, **rules):
    """
    substitution

    >>> from rmpoly import *
    >>> from fractions import Fraction as mpq
    >>> rp,y,x = rgens('y,x',8,mpq)
    >>> p = x**2 + y**2
    >>> p.subs(x=x+y,y=x-y)
     +2*x^2 +2*y^2
    """
    rp = p.rp
    sb = Subs(rp,rp,rules)
    return sb.subs(p)

  def subs_trunc(p, iv,nv, **rules):
    """
    substitution with truncation

    >>> from rmpoly import *
    >>> from fractions import Fraction as mpq
    >>> rp,y,x = rgens('y,x',8,mpq)
    >>> p = x**2 + y**2
    >>> p.subs_trunc(['y','x'],[3,5], y=(x+y)**2)
     +x^4 +4*x^3*y +6*x^2*y^2 +x^2
    """
    rp = p.rp
    sb = Subs(rp,rp,rules)
    return sb.subs_trunc(p,iv,nv)

############## elementary functions ####################

  def _nth_root1(p,n,iv,nv):
    """univariate series nth root of p on commuting ring
    n  integer; compute p^(1/n)
    iv name of the series variable
    nv precision of the series

    The Newton method is used.
    """
    rp = p.rp
    if 0 not in p:
      raise ValueError('no constant term in series')
    if n != int(n):
      raise ValueError('n must be an integer')
    else:
       n = int(n)
    if p[0] != 1:
      raise NotImplementedError('the constant term must be 1')
    p1 = rp(1)
    if p == 1:
      return p
    if n == 0:
      if p != 0:
        return rp('1')
      else:
        raise ValueError('0^0 expression')
    if n == 1:
      return p
    if n < 0:
      n = -n
      sign = 1
    else:
      sign = 0
    for prec in newton_method_sizes(nv):
        tmp = p1.pow_trunc(n+1, iv, prec)
        tmp = tmp.mul_trunc(p, iv, prec)
        p1 += p1/n - tmp/n
    if sign:
      return p1
    else:
      return p1._series_inversion1(iv,nv)
    return p1

  def nth_root(p,n,iv,nv):
    """multivariate series nth root of p
    n  integer; compute p^(1/n)
    iv list of variable names or variable indices
    nv list of truncations for these variables
    In the case of one variable it can also be:
      iv variable name or variable index (0)
      nv truncation integer for the variable
    p is a series with O(x_1^n_1*..x_m^n_m) in 
    variables x_k with index or name iv[k-1]
    p has constant term equal to 1
    TODO: case of constant term different from 1
    >>> from rmpoly import *
    >>> from fractions import Fraction as mpq
    >>> rp,a,x = rgens('a,x',6,mpq)
    >>> (1+a*x).nth_root(3,'x',5)
     -10/243*x^4*a^4 +5/81*x^3*a^3 -1/9*x^2*a^2 +1/3*x*a +1
    """
    rp = p.rp
    # TODO p-1 must not have constant term
    if (p-1).has_constant_term(iv):
      raise NotImplementedError, 'p-1 must not have a constant term in the series variables'
    if hasattr(nv,'__hex__') and rp.commuting:
      return p._nth_root1(n,iv, nv)
    return p.fun('_nth_root1',n,iv,nv)

  def sqrt(p,iv,nv):
    """square root of multivariate series p
    iv list of variable names or variable indices
    nv list of truncations for these variables
    In the case of one variable it can also be:
      iv variable name or variable index (0)
      nv truncation integer for the variable
    p is a series with O(x_1^n_1*..x_m^n_m) in 
    variables x_k with index or name iv[k-1]
    p has constant term equal to 1
    >>> from rmpoly import *
    >>> from fractions import Fraction as mpq
    >>> rp,a,x = rgens('a,x',6,mpq)
    >>> (1+a*x).sqrt('x',5)
     -5/128*x^4*a^4 +1/16*x^3*a^3 -1/8*x^2*a^2 +1/2*x*a +1
    """
    p1 = p.nth_root(-2,iv,nv)
    return p.mul_trunc(p1,iv,nv)



  def _log_series0(p, iv,nv):
    p = p - 1
    rp = p.rp
    s = rp('0')
    q = rp('-1')
    n = 1
    z = rp('0')
    while 1:
      q = -q.mul_trunc(p,iv,nv)
      if q == z:
        break
      s += q/n
      n += 1
    return s

  def _log_series(p, iv,nv):
    if len(p) < 50:
      return p._log_series0(iv,nv)
    rp = p.rp
    if rp.base_field:
      one = rp.base_field(1)
    else:
      one = rp.field(1)
    c = [one*0]
    if hasattr(nv,'__hex__'):
      m = nv
    else:
      m = sum(nv)
    for k in range(1,m):
      if k%2 == 0:
        cf = -one/k
      else:
        cf = one/k
      c.append(cf)
    return (p-1).series_from_list(c,iv,nv)

  def _atan_series0(p, iv,nv):
    rp = p.rp
    s = rp('1')
    q = rp('1')
    n = 3
    p2 = p.square_trunc(iv,nv)
    while 1:
      q = -q.mul_trunc(p2,iv,nv)
      if not q:
        break
      s += q/n
      n += 2
    s = s.mul_trunc(p,iv,nv)
    return s

  def _atan_series(p, iv,nv):
    if len(p) < 50:
      return p._atan_series0(iv,nv)
    rp = p.rp
    mo=p.rp.field(-1)
    c = [-mo]
    if hasattr(nv,'__hex__'):
      m = nv
    else:
      m = sum(nv)
    p2 = p.square_trunc(iv,nv)
    for k in range(1,m):
      c.append(mo**k/(2*k+1))
    s = p2.series_from_list(c,iv,nv)
    s = s.mul_trunc(p,iv,nv)
    #s += p.trunc(iv,nv)
    return s
    
  def _atanh_series0(p, iv,nv):
    rp = p.rp
    s = rp(1)
    q = rp(1)
    n = 3
    p2 = p.square_trunc(iv,nv)
    while 1:
      q = q.mul_trunc(p2,iv,nv)
      if not q:
        break
      s += q/n
      n += 2
    s = s.mul_trunc(p,iv,nv)
    return s

  def _atanh_series(p, iv,nv):
    if len(p) < 50:
      return p._atanh_series0(iv,nv)
    rp = p.rp
    one = rp.field(1)
    c = [one]
    if hasattr(nv,'__hex__'):
      m = nv
    else:
      m = sum(nv)
    p2 = p.square_trunc(iv,nv)
    for k in range(1,m):
      c.append(one/(2*k+1))
    s = p2.series_from_list(c,iv,nv)
    s = s.mul_trunc(p,iv,nv)
    return s

  def log(p,iv,nv):
    """
    logarithm of p with truncation 
    p polynomial starting with 1
    
    For univariate series or with the total degree
    truncation integral dx p^-1*d p/dx is used.
    >>> from rmpoly import *
    >>> from fractions import Fraction as mpq
    >>> rp,y,x = rgens('y,x',10,mpq,order='grlex')
    >>> (1 + x + x*y).log('x',3)
     -1/2*x^2*y^2 -x^2*y -1/2*x^2 +x*y +x
    >>> (1 + x + x*y).log('_t',3)
     -1/2*x^2 +x*y +x
    >>> (1 + x + x*y).log(['y','x'],[2,4])
     +x^3*y +1/3*x^3 -x^2*y -1/2*x^2 +x*y +x
    """
    if (p-1).has_constant_term(iv):
      raise NotImplementedError, 'p-1 must not have a constant term in the series variables'
    if iv in p.rp.pol_gens and hasattr(nv,'__hex__') and p.rp.commuting:
      dlog = p.derivative(iv)
      dlog = dlog.mul_trunc(p._series_inversion1(iv,nv),iv,nv-1)
      return dlog.integrate(iv)
    if p.rp.order == 'grlex' and iv == '_t' and p.rp.commuting:
      p1 = p.series_inversion(iv,nv)
      p2 = p.derivative_t(1)
      dlog = p1.mul_trunc(p2,iv,nv)
      return dlog.derivative_t(-1)
    else:
      return p._log_series(iv,nv)

  def atan(p,iv,nv):
    """
    arctangent of a series
    >>> from rmpoly import *
    >>> from fractions import Fraction as mpq
    >>> rp,a,x = rgens('a,x',6,mpq)
    >>> (a*x).atan('x',10)
     +1/9*x^9*a^9 -1/7*x^7*a^7 +1/5*x^5*a^5 -1/3*x^3*a^3 +x*a
    """
    if p.has_constant_term(iv):
      raise NotImplementedError, 'polynomial must not have constant term in the series variables'
    if iv in p.rp.pol_gens and hasattr(nv,'__hex__') and p.rp.commuting:
      dp = p.derivative(iv)
      p1 = p.square_trunc(iv,nv) + 1
      p1 = p1.series_inversion(iv,nv-1)
      p1 = dp.mul_trunc(p1,iv,nv-1)
      return p1.integrate(iv)
    else:
      return p._atan_series(iv,nv)

  def atanh(p,iv,nv):
    if p.has_constant_term(iv):
      raise NotImplementedError, 'polynomial must not have constant term in the series variables'
    if iv in p.rp.pol_gens and hasattr(nv,'__hex__') and p.rp.commuting:
      dp = p.derivative(iv)
      p1 = -p.square_trunc(iv,nv) + 1
      p1 = p1.series_inversion(iv,nv-1)
      p1 = dp.mul_trunc(p1,iv,nv-1)
      return p1.integrate(iv)
    else:
      return p._atanh_series(iv,nv)

  def _tan1(p,iv,nv):
    rp = p.rp
    p1 = rp(0)
    for prec in newton_method_sizes(nv):
      tmp = p - p1.atan(iv, prec)
      tmp = tmp.mul_trunc(1 + p1.square(),iv, prec)
      p1 += tmp
    return p1

  def tan(p,iv,nv):
    """tangent of a series
    >>> from rmpoly import *
    >>> from fractions import Fraction as mpq
    >>> rp,a,x = rgens('a,x',6,mpq)
    >>> (a*x).tan('x',10)
     +62/2835*x^9*a^9 +17/315*x^7*a^7 +2/15*x^5*a^5 +1/3*x^3*a^3 +x*a
    """
    rp = p.rp
    if p.has_constant_term(iv):
      raise NotImplementedError, 'p must not have constant part in series variables'
    p1 = rp(0)
    if hasattr(nv,'__hex__') and rp.commuting:
      if iv in rp.pol_gens and hasattr(nv,'__hex__') or \
        rp.order == 'grlex' and iv == '_t':
        return p._tan1(iv, nv)
    if hasattr(nv,'__hex__'):
      nv = [int(nv)]
    return p.fun('tan',iv,nv)


  def _tanh1(p,iv,nv):
    rp = p.rp
    p1 = rp('0')
    for prec in newton_method_sizes(nv):
      tmp = p - p1.atanh(iv, prec)
      tmp = tmp.mul_trunc(1 - p1.square(),iv, prec)
      p1 += tmp
    return p1


  def tanh(p,iv,nv):
    """ hyperbolic tangent of a series  
    >>> from rmpoly import *
    >>> from fractions import Fraction as mpq
    >>> rp,a,x = rgens('a,x',6,mpq)
    >>> (a*x).tanh('x',8)
     -17/315*x^7*a^7 +2/15*x^5*a^5 -1/3*x^3*a^3 +x*a
    """
    rp = p.rp
    if p.has_constant_term(iv):
      raise NotImplementedError, 'p must not have constant part in series variables'
    p1 = rp('0')
    if hasattr(nv,'__hex__') and rp.commuting:
      if iv in rp.pol_gens  or \
          rp.order == 'grlex' and iv == '_t':
        return p._tanh1(iv, nv)
    if hasattr(nv,'__hex__'):
      nv = [int(nv)]
    return p.fun('tanh',iv,nv)


    
  def _exp_series0(p,iv,nv):
    if 0 in p:
      raise NotImplementedError, 'p must not have constant part'
    rp = p.rp
    q = rp('1')
    n = 1
    k = 1
    p1 = rp('1')
    z = rp('0')
    while 1:
      q = q.mul_trunc(p,iv,nv)
      if not q:
        break
      p1 += q/n
      k += 1
      n *= k
    return p1


  def _exp_series(p,iv,nv):
    if 0 in p:
      raise NotImplementedError, 'p must not have constant part'
    n=p.rp.field(1)
    k=1
    c = []
    if hasattr(nv,'__hex__'):
      m = nv
    else:
      m = sum(nv)
    for k in range(m):
      c.append(1/n)
      k += 1
      n *= k
    return p.series_from_list(c,iv,nv)


  def _exp1(p,iv,nv):
    """
    exponential of a univariate series, or of a multivariate
    series with total degree truncation
    """
    rp = p.rp
    # TODO p must not have terms not containing series variables
    if 0 in p:
      raise NotImplementedError, 'p must not have constant part'
    p1 = rp('1')
    for prec in newton_method_sizes(nv):
      tmp = (p - p1.log(iv, prec)).mul_trunc(p1,iv, prec)
      p1 += tmp
    return p1

  def exp(p,iv,nv):
    """
    exponential of a series
    >>> from rmpoly import *
    >>> from fractions import Fraction as mpq
    >>> rp,a,x = rgens('a,x',6,mpq)
    >>> (a*x).exp('x',6)
     +1/120*x^5*a^5 +1/24*x^4*a^4 +1/6*x^3*a^3 +1/2*x^2*a^2 +x*a +1
    """
    rp = p.rp
    if p.has_constant_term(iv):
      raise NotImplementedError, 'p must not have constant part in series variables'
    p1 = rp('1')
    if hasattr(nv,'__hex__') and len(p) > 20 and rp.commuting:
      if iv in rp.pol_gens and hasattr(nv,'__hex__') or \
          rp.order == 'grlex' and iv == '_t':
        return p._exp1(iv, nv)
    if rp.base_field:
      one = rp.base_field(1)
    else:
      one = rp.field(1)
    n=1
    k=1
    c = []
    if hasattr(nv,'__hex__'):
      m = nv
    else:
      m = sum(nv)
    for k in range(m):
      c.append(one/n)
      k += 1
      n *= k
    return p.series_from_list(c,iv,nv)
  
  def sin(p,iv,nv):
    """ sin of a series
    >>> from rmpoly import *
    >>> from fractions import Fraction as mpq
    >>> rp,a,x = rgens('a,x',6,mpq)
    >>> (a*x).sin('x',6)
     +1/120*x^5*a^5 -1/6*x^3*a^3 +x*a
    """
    t = (p/2).tan(iv,nv)
    t2 = t.square_trunc(iv,nv)
    p1 = (1+t2).series_inversion(iv,nv)
    return p1.mul_trunc(2*t,iv,nv)

  def cos(p,iv,nv):
    """ cos of a series
    >>> from rmpoly import *
    >>> from fractions import Fraction as mpq
    >>> rp,a,x = rgens('a,x',6,mpq)
    >>> (a*x).cos('x',6)
     +1/24*x^4*a^4 -1/2*x^2*a^2 +1
    """
    t = (p/2).tan(iv,nv)
    t2 = t.square_trunc(iv,nv)
    p1 = (1+t2).series_inversion(iv,nv)
    return p1.mul_trunc(1-t2,iv,nv)

  def cos_sin(p,iv,nv):
    t = (p/2).tan(iv,nv)
    t2 = t.square_trunc(iv,nv)
    p1 = (1+t2).series_inversion(iv,nv)
    return (p1.mul_trunc(1-t2,iv,nv), p1.mul_trunc(2*t,iv,nv))

  def sinh(p,iv,nv):
    """ hyperbolic sin of a series
    >>> from rmpoly import *
    >>> from fractions import Fraction as mpq
    >>> rp,a,x = rgens('a,x',6,mpq)
    >>> (a*x).sinh('x',6)
     +1/120*x^5*a^5 +1/6*x^3*a^3 +x*a
    """
    t = p.exp(iv,nv)
    t1 = t.series_inversion(iv,nv)
    return (t - t1)/2

  def cosh(p,iv,nv):
    """ cos of a series
    >>> from rmpoly import *
    >>> from fractions import Fraction as mpq
    >>> rp,a,x = rgens('a,x',6,mpq)
    >>> (a*x).cosh('x',6)
     +1/24*x^4*a^4 +1/2*x^2*a^2 +1
    """
    t = p.exp(iv,nv)
    t1 = t.series_inversion(iv,nv)
    return (t + t1)/2

  def cosh_sinh(p,iv,nv):
    t = p.exp(iv,nv)
    t1 = t.series_inversion(iv,nv)
    return (t + t1)/2, (t - t1)/2


  def exponents(p):
    decode = p.rp.decode_expv
    ngens = p.rp.ngens
    exps = [decode(expv) for expv in p]
    exps = []
    for expv in p:
      exp = decode(expv)
      n = len(exp)
      if n < ngens:
        exps.append(exp+[0]*(ngens-n))
      else:
        exps.append(exp)
    return exps

  def coefficients(p):
    return p.values()


class Subs:
  """class for substitutions of variables with polynomials,
  possibly truncated.
  Within each substitution the powers of the variable substituted
  are cached.

  Substitutions between different polynomial rings is implemented
  only in the case of rational numbers Fraction, mpq converted to mpf
  unless there is some builtin conversion, like from QQ to
  RealField(prec) in Sage
  TODO: improve these cases

  Substitutions in the noncommutative case are not guaranteed to
  be consistent if the replacing polynomial does not commute
  with all coefficients of terms in which the substitution is
  performed.

  Note
  If all substitutions were with the same truncation
  (or with no truncation) one could cache values at class level.
  e.g. one could do it if Subs had only the method subs; then
  multiple substitutions would not require recomputing cached values.
  """
  def __init__(self,rp1,rp2,rules):
    self.rp1 = rp1
    gens_dict = rp1.gens_dict
    self.rp2 = rp2
    self.with_rp1a = 0
    # if rp2.field is different from rp1.field
    if rp1.field != rp2.field:
      # there is no automatic convertion from mpq or Fraction
      # to mpmath.mpf
      sfld1 = str(rp1.field)
      sfld2 = str(rp2.field)
      if 'mpf' in sfld2:
        if 'mpq' in sfld1 or 'Fraction' in sfld1:
          rp1a = RPoly(rp1.pol_gens,rp1.bits_exp,rp2.field,order=rp1.order)
          self.rp1a = rp1a
          self.with_rp1a = 1
    d = {} # replace monomials with (i,pw)
    if rp1 == rp2:
      gens = rp1.gens()
      for i in range(rp1.ngens):
        d[(i,1)] = gens[i]
    for var in rules:
      d[(gens_dict[var],1)] = rules[var]
    self.d = d

  def q2mpf(self,p):
    """
    
    """
    rp1 = self.rp1
    rp2 = self.rp2
    mpf = rp2.field
    p1 = self.rp1a(0)
    # in Sage QQ numerator and denominator are not properties
    if type(rp1.field(0).numerator) == type(len):
      for expv, val in p.iteritems():
        v = mpf(val.numerator())
        v /= val.denominator()
        p1[expv] = v
    # assume that the rational number has properties
    # numerator and denominator, as in Fraction and
    # in gmpy with version >= 1.14
    else:
      for expv, val in p.iteritems():
        v = mpf(val.numerator)
        v /= val.denominator
        p1[expv] = v
    return p1

  def subs(self,p):
    rp1 = self.rp1
    rp2 = self.rp2
    assert p.rp == self.rp1
    if self.with_rp1a:
      p = self.q2mpf(p)
      rp1 = self.rp1a
    d = self.d.copy()
    p1 = Poly(rp2)
    for expv in p:
      p2 = rp2(1)
      pws = rp1.decode_expv(expv)
      for i in range(len(pws)):
        pw = pws[i]
        if pw == 0:
          continue
        if (i,pw) not in d:
          if pw%2 == 0 and (i,pw/2) in d:
            d[(i,pw)] = d[(i,pw/2)]**2
          elif (i,pw-1) in d:
              d[(i,pw)] = d[(i,pw-1)]*d[(i,1)]
          else:
            d[(i,pw)] = d[(i,1)]**pw
        p2 *= d[(i,pw)]
      p1 += p2*p[expv]
    return p1

  def subs_trunc(self,p,ii,h):
    """substitution with truncation of variable(s) corresponding
    to ii and truncation order(s) h
    """
    rp1 = self.rp1
    rp2 = self.rp2
    assert p.rp == self.rp1
    if self.with_rp1a:
      p = self.q2mpf(p)
      rp1 = self.rp1a
    d = self.d.copy()
    p1 = rp2(0)
    for expv in p:
      p2 = rp2(1)
      pws = rp1.decode_expv(expv)
      for i in range(len(pws)):
        pw = pws[i]
        if pw == 0:
          continue
        if (i,pw) not in d:
          if pw%2 == 0 and (i,pw/2) in d:
            d[(i,pw)] = d[(i,pw/2)].square_trunc(ii,h)
          elif (i,pw-1) in d:
              d[(i,pw)] = d[(i,pw-1)].mul_trunc(d[(i,1)],ii,h)
          else:
            d[(i,pw)] = d[(i,1)].pow_trunc(pw,ii,h)
        p2 = p2.mul_trunc(d[(i,pw)],ii,h)
      p1 += p2*p[expv]
    return p1


if __name__ == "__main__":
    import doctest
    import sys
    if sys.version_info < (2, 6):
      print 'doctests require Fraction, available from Python2.6'
      sys.exit()
    doctest.testmod()

